Data analysis for:

Padilla Perez et al., (2022). The correlated evolution of foraging mode and reproductive output in lizards. Proceedings of the Royal Society B.

Dylan J. Padilla Perez, School of Life Sciences, Arizona State University, Tempe, AZ 85287, USA.

Library

require(AICcmodavg)
require(ape)
require(caper)
require(geiger)
require(kableExtra)
require(MuMIn)
require(nlme)
require(phytools)
require(plotrix)
require(shape)
require(car)
R.version
                 _                           
  platform       x86_64-pc-linux-gnu         
  arch           x86_64                      
  os             linux-gnu                   
  system         x86_64, linux-gnu           
  status                                     
  major          4                           
  minor          1.0                         
  year           2021                        
  month          05                          
  day            18                          
  svn rev        80317                       
  language       R                           
  version.string R version 4.1.0 (2021-05-18)
  nickname       Camp Pontanezen
sessionInfo()
  R version 4.1.0 (2021-05-18)
  Platform: x86_64-pc-linux-gnu (64-bit)
  Running under: CentOS Linux 7 (Core)
  
  Matrix products: default
  BLAS:   /packages/7x/R/R-4.1.0/lib/libRblas.so
  LAPACK: /packages/7x/R/R-4.1.0/lib/libRlapack.so
  
  locale:
   [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
   [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
   [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
   [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
   [9] LC_ADDRESS=C               LC_TELEPHONE=C            
  [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       
  
  attached base packages:
  [1] stats     graphics  grDevices utils     datasets  methods   base     
  
  other attached packages:
   [1] car_3.0-13       carData_3.0-5    shape_1.4.6      plotrix_3.8-2   
   [5] phytools_1.0-1   maps_3.4.0       nlme_3.1-153     MuMIn_1.43.17   
   [9] kableExtra_1.3.4 geiger_2.0.7     caper_1.0.1      mvtnorm_1.1-3   
  [13] MASS_7.3-55      ape_5.6-2        AICcmodavg_2.3-1 knitr_1.37      
  [17] rmarkdown_2.13  
  
  loaded via a namespace (and not attached):
   [1] Rcpp_1.0.8.3            subplex_1.7             svglite_2.1.0          
   [4] lattice_0.20-45         digest_0.6.29           unmarked_1.1.1         
   [7] R6_2.5.1                plyr_1.8.7              stats4_4.1.0           
  [10] evaluate_0.15           coda_0.19-4             httr_1.4.2             
  [13] rlang_0.4.12            rstudioapi_0.13         raster_3.5-15          
  [16] jquerylib_0.1.4         phangorn_2.8.1          Matrix_1.4-0           
  [19] combinat_0.0-8          splines_4.1.0           webshot_0.5.3          
  [22] stringr_1.4.0           munsell_0.5.0           igraph_1.2.11          
  [25] compiler_4.1.0          numDeriv_2016.8-1.1     xfun_0.30              
  [28] systemfonts_1.0.4       pkgconfig_2.0.3         mnormt_2.0.2           
  [31] tmvnsim_1.0-2           htmltools_0.5.2         expm_0.999-6           
  [34] quadprog_1.5-8          codetools_0.2-18        viridisLite_0.4.0      
  [37] grid_4.1.0              lifecycle_1.0.1         jsonlite_1.8.0         
  [40] xtable_1.8-4            magrittr_2.0.3          scales_1.1.1           
  [43] stringi_1.7.6           scatterplot3d_0.3-41    sp_1.4-6               
  [46] xml2_1.3.3              bslib_0.3.1             fastmatch_1.1-3        
  [49] deSolve_1.30            tools_4.1.0             glue_1.6.2             
  [52] abind_1.4-5             parallel_4.1.0          fastmap_1.1.0          
  [55] survival_3.3-1          yaml_2.3.5              colorspace_2.0-3       
  [58] terra_1.5-21            rvest_1.0.2             VGAM_1.1-5             
  [61] clusterGeneration_1.3.7 sass_0.4.1
evo_bd <- read.csv("evo_body_size_mass2.csv", row.names = 1)
evo_bd <- evo_bd[evo_bd$diet != "Herbivorous", ]
evo_bd <- evo_bd[evo_bd$Family != "Sphenodontidae", ]
tree <- read.nexus("time.tree.nex")
str(evo_bd)
  'data.frame': 588 obs. of  25 variables:
   $ reproductive.mode                         : chr  "Oviparous" "Oviparous" "Oviparous" "Oviparous" ...
   $ Genus                                     : chr  "Acanthocercus" "Acanthodactylus" "Acanthodactylus" "Acanthodactylus" ...
   $ Family                                    : chr  "Agamidae" "Lacertidae" "Lacertidae" "Lacertidae" ...
   $ female.SVL                                : num  113.9 64.9 60.1 68.6 61 ...
   $ hatchling.neonate.SVL                     : num  34 25 26.5 29 30 35.8 24.3 65.3 47 27 ...
   $ female.log                                : num  2.1 1.8 1.8 1.8 1.8 1.8 1.7 2.1 2 1.9 ...
   $ hatch.log                                 : num  1.5 1.4 1.4 1.5 1.5 1.6 1.4 1.8 1.7 1.4 ...
   $ intercept                                 : num  -4.69 -4.54 -4.54 -4.54 -4.54 ...
   $ slope                                     : num  3.1 2.95 2.95 2.95 2.95 ...
   $ female_mass                               : num  349 187 173 198 176 ...
   $ hatchling_mass                            : num  100.9 69.2 73.7 81 84 ...
   $ femass.log                                : num  1.7 0.8 0.7 0.9 0.7 0.9 0.5 1.6 1.4 1.3 ...
   $ hatcmass.log                              : num  0.07 -0.42 -0.34 -0.23 -0.18 0.04 -0.45 0.74 0.24 -0.24 ...
   $ rawFem                                    : num  50.06 6.38 5.09 7.52 5.32 ...
   $ rawHatc                                   : num  1.17 0.38 0.45 0.59 0.65 1.1 0.35 5.44 1.74 0.57 ...
   $ largest.clutch                            : int  17 8 8 8 7 8 6 1 2 18 ...
   $ smallest.clutch                           : int  1 3 1 1 2 1 1 1 1 8 ...
   $ largest.mean.clutch.size                  : num  11.4 5 6 4.4 4.8 3 3.3 1 2 14 ...
   $ smallest.mean.clutch.size                 : num  11.4 4.8 2.6 3 2.5 2.6 2.5 1 2 11.5 ...
   $ Longitude.centroid..from.Roll.et.al..2017.: num  30.3 34.8 19.8 -3.2 25.1 ...
   $ Latitude.centroid..from.Roll.et.al..2017. : num  -9.53 31.17 23.35 37.24 31.05 ...
   $ Activity.time                             : chr  "Diurnal" "Diurnal" "Diurnal" "Diurnal" ...
   $ foraging.mode                             : chr  "sit and wait" "widely foraging" "widely foraging" "mixed" ...
   $ diet                                      : chr  "Carnivorous" "Carnivorous" "Carnivorous" "Carnivorous" ...
   $ Binomial                                  : chr  "Acanthocercus_atricollis" "Acanthodactylus_beershebensis" "Acanthodactylus_boskianus" "Acanthodactylus_erythrurus" ...
## Scaled mass index for females

par(mar = c(6, 6, 1, 1), mgp = c(4, 1, 0))
plot(rawFem ~ female.SVL, data = evo_bd, ylab = "Female mass (g)", xlab = "Female length (mm)", las = 1)

par(mar = c(5.1, 4.1, 4.1, 2.1), mgp = c(3, 1, 0))
plot(log(rawFem) ~ log(female.SVL), data = evo_bd, ylab = "Female mass ln(g)", xlab = "Female length ln(mm)", las = 1)

step1 <- lm(log(rawFem) ~ log(female.SVL), data = evo_bd)
summary(step1)
  
  Call:
  lm(formula = log(rawFem) ~ log(female.SVL), data = evo_bd)
  
  Residuals:
      Min      1Q  Median      3Q     Max 
  -2.6266 -0.1464 -0.0359  0.2153  0.7334 
  
  Coefficients:
                   Estimate Std. Error t value Pr(>|t|)    
  (Intercept)     -10.51927    0.11602  -90.67   <2e-16 ***
  log(female.SVL)   2.96706    0.02686  110.48   <2e-16 ***
  ---
  Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
  
  Residual standard error: 0.3334 on 586 degrees of freedom
  Multiple R-squared:  0.9542,  Adjusted R-squared:  0.9541 
  F-statistic: 1.221e+04 on 1 and 586 DF,  p-value: < 2.2e-16
L0 <- mean(evo_bd$female.SVL)
Mi <- evo_bd$rawFem
Li <- evo_bd$female.SVL

evo_bd$M_females <- Mi*(L0/Li)^3


## Scaled mass index for hatchlings

plot(rawHatc ~ hatchling.neonate.SVL, data = evo_bd, ylab = "Hatchling mass (g)", xlab = "Hatchling length (mm)", las = 1)

plot(log(rawHatc) ~ log(hatchling.neonate.SVL), data = evo_bd, ylab = "Hatchling mass log(g)", xlab = "Hatchling length log(mm)", las = 1)

step1h <- lm(log(rawHatc) ~ log(hatchling.neonate.SVL), data = evo_bd)
summary(step1h)
  
  Call:
  lm(formula = log(rawHatc) ~ log(hatchling.neonate.SVL), data = evo_bd)
  
  Residuals:
       Min       1Q   Median       3Q      Max 
  -1.70359 -0.30272  0.02864  0.17771  2.09892 
  
  Coefficients:
                              Estimate Std. Error t value Pr(>|t|)    
  (Intercept)                -10.11741    0.14347  -70.52   <2e-16 ***
  log(hatchling.neonate.SVL)   2.82962    0.04232   66.87   <2e-16 ***
  ---
  Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
  
  Residual standard error: 0.4173 on 586 degrees of freedom
  Multiple R-squared:  0.8841,  Adjusted R-squared:  0.8839 
  F-statistic:  4471 on 1 and 586 DF,  p-value: < 2.2e-16
L0h <- mean(evo_bd$hatchling.neonate.SVL)
Mih <- evo_bd$rawHatc
Lih <- evo_bd$hatchling.neonate.SVL

evo_bd$M_hatchlings <- Mih*(L0h/Lih)^2.8
head(evo_bd)
                                reproductive.mode           Genus     Family
  Acanthocercus_atricollis              Oviparous   Acanthocercus   Agamidae
  Acanthodactylus_beershebensis         Oviparous Acanthodactylus Lacertidae
  Acanthodactylus_boskianus             Oviparous Acanthodactylus Lacertidae
  Acanthodactylus_erythrurus            Oviparous Acanthodactylus Lacertidae
  Acanthodactylus_pardalis              Oviparous Acanthodactylus Lacertidae
  Acanthodactylus_schreiberi            Oviparous Acanthodactylus Lacertidae
                                female.SVL hatchling.neonate.SVL female.log
  Acanthocercus_atricollis           113.9                  34.0        2.1
  Acanthodactylus_beershebensis       64.9                  25.0        1.8
  Acanthodactylus_boskianus           60.1                  26.5        1.8
  Acanthodactylus_erythrurus          68.6                  29.0        1.8
  Acanthodactylus_pardalis            61.0                  30.0        1.8
  Acanthodactylus_schreiberi          69.4                  35.8        1.8
                                hatch.log intercept slope female_mass
  Acanthocercus_atricollis            1.5    -4.686 3.105       349.0
  Acanthodactylus_beershebensis       1.4    -4.543 2.951       187.0
  Acanthodactylus_boskianus           1.4    -4.543 2.951       172.8
  Acanthodactylus_erythrurus          1.5    -4.543 2.951       197.9
  Acanthodactylus_pardalis            1.5    -4.543 2.951       175.5
  Acanthodactylus_schreiberi          1.6    -4.543 2.951       200.3
                                hatchling_mass femass.log hatcmass.log rawFem
  Acanthocercus_atricollis               100.9        1.7         0.07  50.06
  Acanthodactylus_beershebensis           69.2        0.8        -0.42   6.38
  Acanthodactylus_boskianus               73.7        0.7        -0.34   5.09
  Acanthodactylus_erythrurus              81.0        0.9        -0.23   7.52
  Acanthodactylus_pardalis                84.0        0.7        -0.18   5.32
  Acanthodactylus_schreiberi             101.1        0.9         0.04   7.78
                                rawHatc largest.clutch smallest.clutch
  Acanthocercus_atricollis         1.17             17               1
  Acanthodactylus_beershebensis    0.38              8               3
  Acanthodactylus_boskianus        0.45              8               1
  Acanthodactylus_erythrurus       0.59              8               1
  Acanthodactylus_pardalis         0.65              7               2
  Acanthodactylus_schreiberi       1.10              8               1
                                largest.mean.clutch.size
  Acanthocercus_atricollis                          11.4
  Acanthodactylus_beershebensis                      5.0
  Acanthodactylus_boskianus                          6.0
  Acanthodactylus_erythrurus                         4.4
  Acanthodactylus_pardalis                           4.8
  Acanthodactylus_schreiberi                         3.0
                                smallest.mean.clutch.size
  Acanthocercus_atricollis                           11.4
  Acanthodactylus_beershebensis                       4.8
  Acanthodactylus_boskianus                           2.6
  Acanthodactylus_erythrurus                          3.0
  Acanthodactylus_pardalis                            2.5
  Acanthodactylus_schreiberi                          2.6
                                Longitude.centroid..from.Roll.et.al..2017.
  Acanthocercus_atricollis                                           30.28
  Acanthodactylus_beershebensis                                      34.83
  Acanthodactylus_boskianus                                          19.84
  Acanthodactylus_erythrurus                                         -3.20
  Acanthodactylus_pardalis                                           25.10
  Acanthodactylus_schreiberi                                         33.23
                                Latitude.centroid..from.Roll.et.al..2017.
  Acanthocercus_atricollis                                          -9.53
  Acanthodactylus_beershebensis                                     31.17
  Acanthodactylus_boskianus                                         23.35
  Acanthodactylus_erythrurus                                        37.24
  Acanthodactylus_pardalis                                          31.05
  Acanthodactylus_schreiberi                                        35.05
                                Activity.time   foraging.mode        diet
  Acanthocercus_atricollis            Diurnal    sit and wait Carnivorous
  Acanthodactylus_beershebensis       Diurnal widely foraging Carnivorous
  Acanthodactylus_boskianus           Diurnal widely foraging Carnivorous
  Acanthodactylus_erythrurus          Diurnal           mixed Carnivorous
  Acanthodactylus_pardalis            Diurnal    sit and wait Carnivorous
  Acanthodactylus_schreiberi          Diurnal widely foraging Carnivorous
                                                     Binomial M_females
  Acanthocercus_atricollis           Acanthocercus_atricollis  22.23411
  Acanthodactylus_beershebensis Acanthodactylus_beershebensis  15.31746
  Acanthodactylus_boskianus         Acanthodactylus_boskianus  15.38844
  Acanthodactylus_erythrurus       Acanthodactylus_erythrurus  15.28782
  Acanthodactylus_pardalis           Acanthodactylus_pardalis  15.38233
  Acanthodactylus_schreiberi       Acanthodactylus_schreiberi  15.27571
                                M_hatchlings
  Acanthocercus_atricollis         0.9884380
  Acanthodactylus_beershebensis    0.7593747
  Acanthodactylus_boskianus        0.7638861
  Acanthodactylus_erythrurus       0.7781121
  Acanthodactylus_pardalis         0.7796118
  Acanthodactylus_schreiberi       0.8043124
check <- name.check(tree, evo_bd)
rm_phy <- check$tree_not_data
rm_dat <- check$data_not_tree
ctree <- drop.tip(tree, rm_phy)

cdat <- subset(evo_bd, subset = evo_bd$Binomial %in% ctree$tip, select = names(evo_bd)[1:27])
name.check(ctree, cdat)
  [1] "OK"
## Mean clutch size

cdat$mclutch <- (cdat$smallest.mean.clutch.size + cdat$largest.mean.clutch.size) / 2


## Reproductive output

cdat$rep_out <- (cdat$mclutch * cdat$M_hatchlings)
str(cdat)
  'data.frame': 485 obs. of  29 variables:
   $ reproductive.mode                         : chr  "Oviparous" "Oviparous" "Oviparous" "Oviparous" ...
   $ Genus                                     : chr  "Acanthocercus" "Acanthodactylus" "Acanthodactylus" "Acanthodactylus" ...
   $ Family                                    : chr  "Agamidae" "Lacertidae" "Lacertidae" "Lacertidae" ...
   $ female.SVL                                : num  113.9 64.9 60.1 68.6 61 ...
   $ hatchling.neonate.SVL                     : num  34 25 26.5 29 30 35.8 24.3 65.3 47 27 ...
   $ female.log                                : num  2.1 1.8 1.8 1.8 1.8 1.8 1.7 2.1 2 1.9 ...
   $ hatch.log                                 : num  1.5 1.4 1.4 1.5 1.5 1.6 1.4 1.8 1.7 1.4 ...
   $ intercept                                 : num  -4.69 -4.54 -4.54 -4.54 -4.54 ...
   $ slope                                     : num  3.1 2.95 2.95 2.95 2.95 ...
   $ female_mass                               : num  349 187 173 198 176 ...
   $ hatchling_mass                            : num  100.9 69.2 73.7 81 84 ...
   $ femass.log                                : num  1.7 0.8 0.7 0.9 0.7 0.9 0.5 1.6 1.4 1.3 ...
   $ hatcmass.log                              : num  0.07 -0.42 -0.34 -0.23 -0.18 0.04 -0.45 0.74 0.24 -0.24 ...
   $ rawFem                                    : num  50.06 6.38 5.09 7.52 5.32 ...
   $ rawHatc                                   : num  1.17 0.38 0.45 0.59 0.65 1.1 0.35 5.44 1.74 0.57 ...
   $ largest.clutch                            : int  17 8 8 8 7 8 6 1 2 18 ...
   $ smallest.clutch                           : int  1 3 1 1 2 1 1 1 1 8 ...
   $ largest.mean.clutch.size                  : num  11.4 5 6 4.4 4.8 3 3.3 1 2 14 ...
   $ smallest.mean.clutch.size                 : num  11.4 4.8 2.6 3 2.5 2.6 2.5 1 2 11.5 ...
   $ Longitude.centroid..from.Roll.et.al..2017.: num  30.3 34.8 19.8 -3.2 25.1 ...
   $ Latitude.centroid..from.Roll.et.al..2017. : num  -9.53 31.17 23.35 37.24 31.05 ...
   $ Activity.time                             : chr  "Diurnal" "Diurnal" "Diurnal" "Diurnal" ...
   $ foraging.mode                             : chr  "sit and wait" "widely foraging" "widely foraging" "mixed" ...
   $ diet                                      : chr  "Carnivorous" "Carnivorous" "Carnivorous" "Carnivorous" ...
   $ Binomial                                  : chr  "Acanthocercus_atricollis" "Acanthodactylus_beershebensis" "Acanthodactylus_boskianus" "Acanthodactylus_erythrurus" ...
   $ M_females                                 : num  22.2 15.3 15.4 15.3 15.4 ...
   $ M_hatchlings                              : num  0.988 0.759 0.764 0.778 0.78 ...
   $ mclutch                                   : num  11.4 4.9 4.3 3.7 3.65 ...
   $ rep_out                                   : num  11.27 3.72 3.28 2.88 2.85 ...
tapply(cdat$M_females, cdat$Family, length)
           Agamidae          Anguidae       Anniellidae  Carphodactylidae 
                 39                 8                 2                 4 
     Chamaeleonidae        Cordylidae    Corytophanidae     Crotaphytidae 
                  5                 6                 3                 5 
        Dactyloidae   Diplodactylidae     Eublepharidae        Gekkonidae 
                 40                14                 6                35 
     Gerrhosauridae  Gymnophthalmidae    Helodermatidae     Hoplocercidae 
                  1                17                 2                 1 
         Lacertidae    Leiocephalidae       Liolaemidae         Opluridae 
                 49                 3                11                 1 
    Phrynosomatidae  Phyllodactylidae     Polychrotidae       Pygopodidae 
                 39                15                 2                 4 
          Scincidae     Shinisauridae Sphaerodactylidae           Teiidae 
                 94                 1                13                20 
       Tropiduridae         Varanidae       Xantusiidae      Xenosauridae 
                 12                28                 3                 2



par(mar = c(1, 1, 1, 1))
plot(NA, NA, xlim = c(0, 10), ylim = c(-1, 10), ann = FALSE, axes = FALSE)


text(5, 10, "Foraging mode")
segments(3.5, 9.6, 6.5, 9.6, lwd = 2)
segments(3.5, 10.4, 6.5, 10.4, lwd = 2)
segments(3.5, 10.4, 3.5, 9.6, lwd = 2)
segments(6.5, 10.4, 6.5, 9.6, lwd = 2)
text(5, 7, "Distance traveled")
segments(3.5, 6.5, 6.5, 6.5, lwd = 2)
segments(3.5, 7.5, 6.5, 7.5, lwd = 2)
segments(3.5, 6.5, 3.5, 7.5, lwd = 2)
segments(6.5, 7.5, 6.5, 6.5, lwd = 2)
text(4.97, 3.7, "Energy gain")
segments(3.7, 4.2, 6.2, 4.2, lwd = 2)
segments(3.7, 3.2, 6.2, 3.2, lwd = 2)
segments(3.7, 4.2, 3.7, 3.2, lwd = 2)
segments(6.2, 4.2, 6.2, 3.2, lwd = 2)
text(8.6, 7, "Predation risk")
segments(7.5, 6.5, 9.8, 6.5, lwd = 2)
segments(7.5, 7.5, 9.8, 7.5, lwd = 2)
segments(7.5, 6.5, 7.5, 7.5, lwd = 2)
segments(9.8, 7.5, 9.8, 6.5, lwd = 2)

text(8.64, 3.7, "Reproductive output")
segments(7.2, 4.2, 10.1, 4.2, lwd = 2)
segments(7.2, 3.2, 10.1, 3.2, lwd = 2)
segments(7.2, 4.2, 7.2, 3.2, lwd = 2)
segments(10.1, 3.2, 10.1, 4.2, lwd = 2)

text(1.57, 3.7, "Energy expenditure")
segments(0, 4.2, 3.1, 4.2, lwd = 2)
segments(0, 3.2, 3.1, 3.2, lwd = 2)
segments(0, 4.2, 0, 3.2, lwd = 2)
segments(3.1, 4.2, 3.1, 3.2, lwd = 2)

text(4.9, 0, "Net energy")
segments(3.7, 0.5, 6.2, 0.5, lwd = 2)
segments(3.7, -0.5, 6.2, -0.5, lwd = 2)
segments(3.7, 0.5, 3.7, -0.5, lwd = 2)
segments(6.2, 0.5, 6.2, -0.5, lwd = 2)

Arrows(5, 9.5, 5, 7.8, lwd = 2, arr.type = "triangle")
Arrows(5, 6.3, 5, 4.8, lwd = 2, arr.type = "triangle")
Arrows(5, 4.8, 5, 6.1, lwd = 2, arr.type = "triangle")
Arrows(3.8, 6.3, 2, 4.8, lwd = 2, arr.type = "triangle")
Arrows(8.3, 4.76, 6.3, 6.2, lwd = 2, arr.type = "triangle")
Arrows(8.7, 6.1, 8.7, 4.8, lwd = 2, arr.type = "triangle")
Arrows(8.7, 4.8, 8.7, 6.1, lwd = 2, arr.type = "triangle")
Arrows(6.6, 7, 7, 7, lwd = 2, arr.type = "triangle")
Arrows(7.1, 3.7, 6.5, 3.7, lwd = 2, arr.type = "triangle")
Arrows(5, 3, 5, 1, lwd = 2, arr.type = "triangle")
Arrows(1.8, 3, 3.8, 1, lwd = 2, arr.type = "triangle")
Arrows(6.1, 1, 8.3, 3, lwd = 2, arr.type = "triangle")

Figure 1. A conceptual model depicting putative relationships among foraging behavior, energetics, predation risk, and reproductive output. The predicted relationships were derived from theoretical models of life-history evolution (see text for details).



layout(matrix(c(0, 1, 0,
                0, 1, 0,
                0, 2, 0,
                0, 2, 0,
                0, 3, 0,
                0, 3, 0), nrow = 6, ncol = 3, byrow = TRUE))


## A


bertalanffy <- function(la, k, t){

    la <- as.numeric(la)
    k <- as.numeric(k)
    t <- as.numeric(t)
    la*(1 - 4*exp(-k*t))

}

fem1 <- bertalanffy(la = 4, k = 0.3, t = seq(1, 35, 1))
fem2 <- bertalanffy(la = 4.9, k = 0.25, t = seq(1, 35, 1))

par(mar = c(3, 3, 2, 1), mgp = c(0.8, 0, 0), cex.lab = 1.2)
plot(fem1 - 0.7 ~ seq(1, 35, 1), xlim = c(-2, 35), ylim = c(-6, 8), type = "l", lwd = 2, xaxt = "n", yaxt = "n", ylab = "Cumulative energy gain (m)", xlab = "time (t)")
lines(seq(1, 35, 1) - 0.48, fem2, lwd = 2, lty = 4)
legend("topleft", legend = c("large female", "small female"), lty = c(4, 1), lwd = c(2, 2), bty = "n", cex = 0.7)
mtext("(a)", at = -6, line = 1, cex = 0.8, font = 2)



## B

fem3 <- bertalanffy(la = 5, k = 0.082, t = seq(1, 35, 1))
fem4 <- bertalanffy(la = 10, k = 0.21, t = seq(1, 35, 1))


par(mar = c(3, 3, 2, 1), mgp = c(0.8, 0, 0))

plot(fem3 - 0.7 ~ seq(1, 35, 1), xlim = c(4, 35), ylim = c(-6, 5), type = "l", lwd = 2,  col = "skyblue", xaxt = "n", yaxt = "n", ylab = "Cumulative energy gain (m)", xlab = " ")
lines(seq(1, 35, 1), fem4 - 8, lwd = 2, lty = 2, col = "purple")
segments(2.6, -6.5, 25, 6.2, lwd = 2)
segments(2.6, -6.5, 33.2, 6, lwd = 2)
segments(13.2, -0.5, 2.6, -0.5, lwd = 2, lty = 2, col = "purple")
segments(15.9, -1.1, 2.6, -1.1, lwd = 2, col = "skyblue")
segments(13.2, -0.5, 13.2, -7, lwd = 2, lty = 2, col = "purple")
segments(15.9, -1.1, 15.9, -7, lwd = 2, col = "skyblue")
legend("topleft", legend = c("sit and wait", "widely foraging"), lty = c(1, 2), lwd = c(2, 2), col = c("skyblue", "purple"), bty = "n", cex = 0.7)
mtext("(b)", at = 0.6, line = 1, cex = 0.8, font = 2)
mtext(expression("t"[wf]), side = 1, at = 13.2, line = 0.5, cex = 0.8, font = 2)
mtext(expression("t"[sw]), side = 1, at = 17.2, line = 0.5, cex = 0.8, font = 2)
mtext("time (t)", side = 1, at = 25, line = 0.9, cex = 0.8)


## C

par(mar = c(3, 3, 2, 1), mgp = c(0.8, 0, 0))

curve(2 - log(x)^2, ylim = c(-0.5, 3), xlim = c(0.1, 1), lwd = 2, lty = 2, col = "purple", xaxt = "n", yaxt = "n", ylab = "Cumulative energy gain (m)", xlab = "time (t)")
curve(log(x) + 1, add = TRUE, lwd = 2, col = "skyblue")
legend("topleft", legend = c("sit and wait", "widely foraging"), lty = c(1, 2), lwd = c(2, 2), col = c("skyblue", "purple"), bty = "n", cex = 0.7)
mtext("(c)", at = 0.001, line = 1, cex = 0.8, font = 2)

Figure 4. A theoretical model relating the cumulative energy gain for reproduction, m, as a function of time spent foraging, t, and maternal size. (a) Larger females reach their maximum capacity at a higher value of m than small females. (b) Widely-foraging females that are smaller but more efficient foragers may produce a greater reproductive output than larger sit-and-wait females. (c) Widely-foraging females may also produce a greater reproductive output than sit-and-wait females if they are both larger and more efficient foragers. twf and tsw in (b) represent the optimal foraging time of widely-foraging females and sit-and-wait females, respectively.




Fitting models (PGLS)


## Lacertilia

cdat$foraging.mode <- factor(cdat$foraging.mode, levels = c("widely foraging", "sit and wait", "mixed")) ## I am changing the order of the levels here so that the contrast of the coefficients can be made against the level "widely foraging".


rawMass <- cdat

rawMass <- rawMass[rawMass$rawFem < 5000, ] # Excluding a mixed-foraging species (too heavy)

check_rawMass <- name.check(tree, rawMass)
rm_phy_rawMass <- check_rawMass$tree_not_data
rm_dat_rawMass <- check_rawMass$data_not_tree
ctree_rawMass <- drop.tip(tree, rm_phy_rawMass)
name.check(ctree_rawMass, rawMass)
  [1] "OK"
## Reproductive output

rawMass$rep_out_mass <- (rawMass$mclutch * rawMass$rawHatc)


## Number of individuals in each category (foraging mode)

tapply(rawMass$rawFem, rawMass$foraging.mode, length)
  widely foraging    sit and wait           mixed 
              202             222              60
## Mean mass per category

tapply(rawMass$rawFem, rawMass$foraging.mode, mean)
  widely foraging    sit and wait           mixed 
        153.48926        32.35770        21.48317
## Mean reproductive output per category

tapply(rawMass$rep_out_mass, rawMass$foraging.mode, mean)
  widely foraging    sit and wait           mixed 
        28.828361        5.540074        6.718067
mod_mass <- gls(rep_out_mass ~ rawFem*foraging.mode, correlation = corBrownian(phy = ctree_rawMass, form = ~Binomial), data = rawMass, method = "ML")
summary(mod_mass)
  Generalized least squares fit by maximum likelihood
    Model: rep_out_mass ~ rawFem * foraging.mode 
    Data: rawMass 
         AIC      BIC    logLik
    4938.271 4967.545 -2462.135
  
  Correlation Structure: corBrownian
   Formula: ~Binomial 
   Parameter estimate(s):
  numeric(0)
  
  Coefficients:
                                       Value Std.Error   t-value p-value
  (Intercept)                       6.963796 27.550624  0.252764  0.8006
  rawFem                            0.139790  0.005549 25.189779  0.0000
  foraging.modesit and wait         0.809812  5.670050  0.142823  0.8865
  foraging.modemixed               -1.512868  6.179807 -0.244808  0.8067
  rawFem:foraging.modesit and wait -0.081644  0.010811 -7.552212  0.0000
  rawFem:foraging.modemixed         0.026696  0.089214  0.299233  0.7649
  
   Correlation: 
                                   (Intr) rawFem frg.aw frgng. rF:.aw
  rawFem                           -0.027                            
  foraging.modesit and wait        -0.105  0.099                     
  foraging.modemixed               -0.063  0.042  0.426              
  rawFem:foraging.modesit and wait  0.022 -0.308 -0.211 -0.088       
  rawFem:foraging.modemixed        -0.004  0.001  0.106 -0.485  0.030
  
  Standardized residuals:
          Min          Q1         Med          Q3         Max 
  -2.21896992 -0.08372645 -0.07286486 -0.05356706  4.38147263 
  
  Residual standard error: 85.60208 
  Degrees of freedom: 484 total; 478 residual
mod_mass1 <- gls(rep_out_mass ~ rawFem + foraging.mode, correlation = corBrownian(phy = ctree_rawMass, form = ~Binomial), data = rawMass, method = "ML")
summary(mod_mass1)
  Generalized least squares fit by maximum likelihood
    Model: rep_out_mass ~ rawFem + foraging.mode 
    Data: rawMass 
         AIC      BIC    logLik
    4989.078 5009.988 -2489.539
  
  Correlation Structure: corBrownian
   Formula: ~Binomial 
   Parameter estimate(s):
  numeric(0)
  
  Coefficients:
                                Value Std.Error   t-value p-value
  (Intercept)               11.599882 29.087362  0.398795  0.6902
  rawFem                     0.126833  0.005574 22.752688  0.0000
  foraging.modesit and wait -8.569257  5.814017 -1.473896  0.1412
  foraging.modemixed        -4.064938  5.688570 -0.714580  0.4752
  
   Correlation: 
                            (Intr) rawFem frg.aw
  rawFem                    -0.021              
  foraging.modesit and wait -0.102  0.036       
  foraging.modemixed        -0.073  0.024  0.545
  
  Standardized residuals:
          Min          Q1         Med          Q3         Max 
  -2.67186452 -0.11700019 -0.04441283 -0.02384972  4.47349177 
  
  Residual standard error: 90.58862 
  Degrees of freedom: 484 total; 480 residual
mod_mass2 <- gls(rep_out_mass ~ rawFem*foraging.mode, correlation = corPagel(value = 0, phy = ctree_rawMass, form = ~Binomial), data = rawMass, method = "ML")
summary(mod_mass2)
  Generalized least squares fit by maximum likelihood
    Model: rep_out_mass ~ rawFem * foraging.mode 
    Data: rawMass 
         AIC      BIC    logLik
    4761.285 4794.742 -2372.642
  
  Correlation Structure: corPagel
   Formula: ~Binomial 
   Parameter estimate(s):
       lambda 
  -0.08292133 
  
  Coefficients:
                                       Value Std.Error   t-value p-value
  (Intercept)                       0.221955  0.383885   0.57818  0.5634
  rawFem                            0.185596  0.002613  71.03444  0.0000
  foraging.modesit and wait         1.803715  0.185985   9.69817  0.0000
  foraging.modemixed               -5.678696  3.737821  -1.51925  0.1294
  rawFem:foraging.modesit and wait -0.122752  0.011515 -10.65994  0.0000
  rawFem:foraging.modemixed         0.106016  0.112114   0.94561  0.3448
  
   Correlation: 
                                   (Intr) rawFem frg.aw frgng. rF:.aw
  rawFem                           -0.279                            
  foraging.modesit and wait         0.168  0.096                     
  foraging.modemixed               -0.907  0.279 -0.482              
  rawFem:foraging.modesit and wait -0.006 -0.473 -0.560  0.021       
  rawFem:foraging.modemixed         0.042 -0.216  0.618 -0.425  0.022
  
  Standardized residuals:
           Min           Q1          Med           Q3          Max 
  -8.754409279 -0.041686637 -0.001960039  0.068491652  9.544845883 
  
  Residual standard error: 32.6597 
  Degrees of freedom: 484 total; 478 residual
Anova(mod_mass2, type = "III")
  Analysis of Deviance Table (Type III tests)
  
  Response: rep_out_mass
                       Df     Chisq Pr(>Chisq)    
  (Intercept)           1    0.3343     0.5631    
  rawFem                1 5045.8917     <2e-16 ***
  foraging.mode         2  107.0707     <2e-16 ***
  rawFem:foraging.mode  2  115.0202     <2e-16 ***
  ---
  Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
mod_mass3 <- gls(rep_out_mass ~ rawFem + foraging.mode, correlation = corPagel(value = 0, phy = ctree_rawMass, form = ~Binomial), data = rawMass)
summary(mod_mass3)
  Generalized least squares fit by REML
    Model: rep_out_mass ~ rawFem + foraging.mode 
    Data: rawMass 
         AIC      BIC    logLik
    4871.726 4896.769 -2429.863
  
  Correlation Structure: corPagel
   Formula: ~Binomial 
   Parameter estimate(s):
        lambda 
  -0.001644695 
  
  Coefficients:
                                Value Std.Error  t-value p-value
  (Intercept)                4.896045  2.608520  1.87694  0.0611
  rawFem                     0.155177  0.004527 34.27451  0.0000
  foraging.modesit and wait -4.181294  3.548787 -1.17823  0.2393
  foraging.modemixed        -1.508168  5.441039 -0.27718  0.7818
  
   Correlation: 
                            (Intr) rawFem frg.aw
  rawFem                    -0.263              
  foraging.modesit and wait -0.740  0.151       
  foraging.modemixed        -0.492  0.110  0.348
  
  Standardized residuals:
           Min           Q1          Med           Q3          Max 
  -8.892572381 -0.119008564 -0.045067474 -0.002331544  9.525337598 
  
  Residual standard error: 36.92721 
  Degrees of freedom: 484 total; 480 residual
## mod_mass2 is the most likely model even when the models below are included in the analysis

#mod_mass4 <- gls(rep_out_mass ~ rawFem*foraging.mode, correlation = corPagel(value = 1, phy = ctree_rawMass, form = ~Binomial), data = rawMass, method = "ML")
#summary(mod_mass4)

#mod_mass5 <- gls(rep_out_mass ~ rawFem + foraging.mode, correlation = corPagel(value = 1, phy = ctree_rawMass, form = ~Binomial), data = rawMass)
#summary(mod_mass5)


## Model selection table

output_rawMass <- model.sel(mod_mass, mod_mass1, mod_mass2, mod_mass3)

sel.table_mass <- round(as.data.frame(output_rawMass)[8:12], 3)
names(sel.table_mass)[1] <- "K"
sel.table_mass$Model <- rownames(sel.table_mass)
sel.table_mass$Model <- as.character(c(formula(mod_mass2), formula(mod_mass3), formula(mod_mass), formula(mod_mass1)))
sel.table_mass <- sel.table_mass[ , c(6, 1, 2, 3, 4, 5)]
sel.table_mass
                                            Model K    logLik     AICc   delta
  mod_mass2 rep_out_mass ~ rawFem * foraging.mode 8 -2372.642 4761.588   0.000
  mod_mass3 rep_out_mass ~ rawFem + foraging.mode 6 -2429.863 4871.902 110.314
  mod_mass  rep_out_mass ~ rawFem * foraging.mode 7 -2462.135 4938.506 176.918
  mod_mass1 rep_out_mass ~ rawFem + foraging.mode 5 -2489.539 4989.203 227.615
            weight
  mod_mass2      1
  mod_mass3      0
  mod_mass       0
  mod_mass1      0




Model checking (raw mass)


plot(mod_mass2, resid(., type = "n") ~ fitted(.), main = "Residuals v Fitted Values", abline = c(0, 0))

res <- resid(mod_mass2, type="n")
qqnorm(res, las = 1)
qqline(res)



A visual check of the most likely model indicates that the homogeneity and normality of residuals seem to be affected by influential cases (i.e., outliers). As a result, we standardized hatchling mass and maternal mass to alleviate this issue (see model checking below). In doing so, we recovered a more symetrical distibution of body mass among lizard species, increasing the robustness of our analysis. Influential cases should not be dropped from the model as there is no reasonable justification to do so. Instead, the patterns observed in our study could be resolved or even alter as more or better data become available.



## Standardized mass, M.

dat_full <- cdat

check_full <- name.check(tree, dat_full)
rm_phy_full <- check_full$tree_not_data
rm_dat_full <- check_full$data_not_tree
ctree_full <- drop.tip(tree, rm_phy_full)
name.check(ctree_full, dat_full)
  [1] "OK"
## Number of individuals in each category (foraging mode)

tapply(dat_full$M_females, dat_full$foraging.mode, length)
  widely foraging    sit and wait           mixed 
              202             222              61
## Mean scaled mass index per category

tapply(dat_full$M_females, dat_full$foraging.mode, mean)
  widely foraging    sit and wait           mixed 
         14.24854        17.88123        15.77951
## Mean reproductive output per category

tapply(dat_full$rep_out, dat_full$foraging.mode, mean)
  widely foraging    sit and wait           mixed 
         3.220662        4.861416        5.788185
full_mod <- gls(rep_out ~ M_females*foraging.mode, correlation = corBrownian(phy = ctree_full, form = ~Binomial), data = dat_full, method = "ML")
summary(full_mod)
  Generalized least squares fit by maximum likelihood
    Model: rep_out ~ M_females * foraging.mode 
    Data: dat_full 
        AIC      BIC   logLik
    2883.82 2913.109 -1434.91
  
  Correlation Structure: corBrownian
   Formula: ~Binomial 
   Parameter estimate(s):
  numeric(0)
  
  Coefficients:
                                          Value Std.Error    t-value p-value
  (Intercept)                          7.394834  3.950072  1.8720755  0.0618
  M_females                           -0.222565  0.143833 -1.5473860  0.1224
  foraging.modesit and wait            2.451992  1.976688  1.2404548  0.2154
  foraging.modemixed                  -1.733346  2.026141 -0.8554912  0.3927
  M_females:foraging.modesit and wait -0.197860  0.121228 -1.6321394  0.1033
  M_females:foraging.modemixed         0.190587  0.126198  1.5102234  0.1316
  
   Correlation: 
                                      (Intr) M_fmls frg.aw frgng. M_:.aw
  M_females                           -0.557                            
  foraging.modesit and wait           -0.193  0.323                     
  foraging.modemixed                  -0.141  0.266  0.442              
  M_females:foraging.modesit and wait  0.170 -0.339 -0.938 -0.500       
  M_females:foraging.modemixed         0.145 -0.309 -0.468 -0.947  0.587
  
  Standardized residuals:
         Min         Q1        Med         Q3        Max 
  -0.8510203 -0.3170531 -0.1719674  0.2257898  5.2981955 
  
  Residual standard error: 10.19489 
  Degrees of freedom: 485 total; 479 residual
full_mod1 <- gls(rep_out ~ M_females + foraging.mode, correlation = corBrownian(phy = ctree_full, form = ~Binomial), data = dat_full, method = "ML")
summary(full_mod1)
  Generalized least squares fit by maximum likelihood
    Model: rep_out ~ M_females + foraging.mode 
    Data: dat_full 
         AIC      BIC    logLik
    2891.793 2912.713 -1440.896
  
  Correlation Structure: corBrownian
   Formula: ~Binomial 
   Parameter estimate(s):
  numeric(0)
  
  Coefficients:
                                Value Std.Error   t-value p-value
  (Intercept)                7.826492  3.926231  1.993386  0.0468
  M_females                 -0.242185  0.135280 -1.790247  0.0740
  foraging.modesit and wait -1.189500  0.662400 -1.795743  0.0732
  foraging.modemixed         1.604453  0.640735  2.504084  0.0126
  
   Correlation: 
                            (Intr) M_fmls frg.aw
  M_females                 -0.537              
  foraging.modesit and wait -0.119  0.062       
  foraging.modemixed        -0.036 -0.046  0.539
  
  Standardized residuals:
         Min         Q1        Med         Q3        Max 
  -0.8269391 -0.2711459 -0.1582826  0.2007470  5.4010162 
  
  Residual standard error: 10.3215 
  Degrees of freedom: 485 total; 481 residual
full_mod2 <- gls(rep_out ~ M_females*foraging.mode, correlation = corPagel(value = 0, phy = ctree_full, form = ~Binomial), data = dat_full, method = "ML")
summary(full_mod2)
  Generalized least squares fit by maximum likelihood
    Model: rep_out ~ M_females * foraging.mode 
    Data: dat_full 
         AIC      BIC    logLik
    2770.838 2804.311 -1377.419
  
  Correlation Structure: corPagel
   Formula: ~Binomial 
   Parameter estimate(s):
     lambda 
  0.8119803 
  
  Coefficients:
                                          Value Std.Error    t-value p-value
  (Intercept)                         -0.194013 2.6555475 -0.0730595  0.9418
  M_females                            0.270669 0.1240541  2.1818666  0.0296
  foraging.modesit and wait            3.101662 1.9832886  1.5638984  0.1185
  foraging.modemixed                  -4.069565 2.2677650 -1.7945267  0.0734
  M_females:foraging.modesit and wait -0.231656 0.1257082 -1.8428111  0.0660
  M_females:foraging.modemixed         0.321048 0.1451586  2.2117076  0.0275
  
   Correlation: 
                                      (Intr) M_fmls frg.aw frgng. M_:.aw
  M_females                           -0.701                            
  foraging.modesit and wait           -0.328  0.452                     
  foraging.modemixed                  -0.223  0.325  0.455              
  M_females:foraging.modesit and wait  0.311 -0.485 -0.949 -0.484       
  M_females:foraging.modemixed         0.221 -0.359 -0.467 -0.957  0.545
  
  Standardized residuals:
         Min         Q1        Med         Q3        Max 
  -1.2090603 -0.3604564 -0.2210199  0.1140136  7.4176637 
  
  Residual standard error: 6.419291 
  Degrees of freedom: 485 total; 479 residual
full_mod3 <- gls(rep_out ~ M_females + foraging.mode, correlation = corPagel(value = 0, phy = ctree_full, form = ~Binomial), data = dat_full, method = "ML")
summary(full_mod3)
  Generalized least squares fit by maximum likelihood
    Model: rep_out ~ M_females + foraging.mode 
    Data: dat_full 
         AIC      BIC    logLik
    2784.802 2809.907 -1386.401
  
  Correlation Structure: corPagel
   Formula: ~Binomial 
   Parameter estimate(s):
     lambda 
  0.8188292 
  
  Coefficients:
                                 Value Std.Error    t-value p-value
  (Intercept)                0.7964322 2.5818235  0.3084766  0.7579
  M_females                  0.2077960 0.1098458  1.8917060  0.0591
  foraging.modesit and wait -0.8114631 0.6242388 -1.2999242  0.1942
  foraging.modemixed         1.1099727 0.6600939  1.6815375  0.0933
  
   Correlation: 
                            (Intr) M_fmls frg.aw
  M_females                 -0.657              
  foraging.modesit and wait -0.122 -0.006       
  foraging.modemixed        -0.080 -0.010  0.513
  
  Standardized residuals:
          Min          Q1         Med          Q3         Max 
  -0.77190136 -0.34187016 -0.21073950  0.06296427  7.81257815 
  
  Residual standard error: 6.59119 
  Degrees of freedom: 485 total; 481 residual
#full_mod4 <- gls(rep_out ~ M_females*foraging.mode, correlation = corPagel(value = 1, phy = ctree_full, form = ~Binomial, fixed = FALSE), data = dat_full, method = "ML")
#summary(full_mod4)

#full_mod5 <- gls(rep_out ~ M_females + foraging.mode, correlation = corPagel(value = 1, phy = ctree_full, form = ~Binomial, fixed = FALSE), data = dat_full, method = "ML")
#summary(full_mod5)

#full_mod6 <- gls(rep_out ~ M_females*foraging.mode, correlation = corMartins(value = 1, phy = ctree_full, form = ~Binomial, fixed = FALSE), data = dat_full, method = "ML")
#summary(full_mod6)




Model checking (scaled mass index)


plot(full_mod2, resid(., type = "n") ~ fitted(.), main = "Residuals v Fitted Values", abline = c(0, 0))

res2 <- resid(full_mod2, type="n")
qqnorm(res2, las = 1)
qqline(res2)




Table S1. Model selection table of candidate models included in the analysis

output_full <- model.sel(full_mod, full_mod1, full_mod2, full_mod3)

sel.table_full <- round(as.data.frame(output_full)[7:11], 3)
names(sel.table_full)[1] <- "K"
sel.table_full$Model <- rownames(sel.table_full)
sel.table_full$Model <- as.character(c(formula(full_mod2), formula(full_mod3), formula(full_mod), formula(full_mod1)))
sel.table_full <- sel.table_full[ , c(6, 1, 2, 3, 4, 5)]


kable(sel.table_full, digits = 3, row.names = FALSE) %>%
  kable_styling(bootstrap_options = "striped", full_width = FALSE, position = "center")
Model K logLik AICc delta weight
rep_out ~ M_females * foraging.mode 8 -1377.419 2771.141 0.000 0.999
rep_out ~ M_females + foraging.mode 6 -1386.401 2784.977 13.837 0.001
rep_out ~ M_females * foraging.mode 7 -1434.910 2884.055 112.915 0.000
rep_out ~ M_females + foraging.mode 5 -1440.896 2891.918 120.777 0.000


Table 1. Analysis of Deviance Table (Type III tests) for the most likely model, based on the ranking of AICc for potential candidate models.

Df Chisq Pr(>Chisq)
(Intercept) 1 0.005 0.942
M_females 1 4.761 0.029
foraging.mode 2 10.372 0.006
M_females:foraging.mode 2 18.105 0.000
layout(matrix(c(0, 0, 0, 0,
                1, 1, 2, 2,
                1, 1, 2, 2,
                3, 3, 4, 4,
                3, 3, 4, 4,
                0, 0, 0, 0), nrow = 6, ncol = 4, byrow = TRUE))


par(mar = c(0, 5, 0, 0.5))


## Distribution of body size


rm1 <- density(rawMass$rawFem[rawMass$foraging.mode == "widely foraging"])
rm2 <- density(rawMass$rawFem[rawMass$foraging.mode == "mixed"])
rm3 <- density(rawMass$rawFem[rawMass$foraging.mode == "sit and wait"])
plot(1, type = "n", ann = FALSE, bty = "n", xaxt = "n", yaxt = "n")


par(mar = c(0, 4, 0, 0.5))

t1 <- density(dat_full$M_females[dat_full$foraging.mode == "widely foraging"])
t2 <- density(dat_full$M_females[dat_full$foraging.mode == "mixed"])
t3 <- density(dat_full$M_females[dat_full$foraging.mode == "sit and wait"])
plot(t1, ylim = c(0, max(c(t1$y, t1$y))), ann = FALSE, bty = "n", xaxt = "n", yaxt = "n")
polygon(t1, col = rgb(0.8, 0, 0.7, alpha = 0.3))
lines(t2)
polygon(t2, col = rgb(0.8, 0, 0, alpha = 0.3))
lines(t3)
polygon(t3, col = rgb(0, 0, 1, alpha = 0.2))



## Lacertilia

## Raw mass


## widely foraging

rw_fg <- rawMass[rawMass$foraging.mode == "widely foraging", ]

check_rm_wf <- name.check(tree, rw_fg)
rm_phy_rm_wf <- check_rm_wf$tree_not_data
rm_dat_rm_wf <- check_rm_wf$data_not_tree
ctree_rm_wf <- drop.tip(tree, rm_phy_rm_wf)
#name.check(ctree_rm_wf, rw_fg)



SSX.rm_wf <- sum(round((rw_fg$rawFem - mean(rw_fg$rawFem))^2), 2)
s2.rm_wf <- var(rw_fg$rep_out_mass)
n.rm_wf <- length(rw_fg$rep_out_mass)
x.rm_wf <- seq(0.50, 4072.0, 1)
m.x.rm_wf <- mean(round(rw_fg$rawFem,1))
se.rm_wf <- sqrt(s2.rm_wf*((1/n.rm_wf) + (((x.rm_wf - m.x.rm_wf)^2)/SSX.rm_wf)))
is.rm_wf <- qt(0.975, df = n.rm_wf - 2)
ii.rm_wf <- qt(0.025, df = n.rm_wf - 2)
ic.s.rm_wf <- se.rm_wf*is.rm_wf
ic.i.rm_wf <- se.rm_wf*ii.rm_wf
upper.i.rm_wf <- (coef(mod_mass2)[1] + coef(mod_mass2)[2]*x.rm_wf) + ic.s.rm_wf
lower.i.rm_wf <- (coef(mod_mass2)[1] + coef(mod_mass2)[2]*x.rm_wf) + ic.i.rm_wf


# sit and wait


rm_sw <- rawMass[rawMass$foraging.mode == "sit and wait", ]

check_rm_sw <- name.check(tree, rm_sw)
rm_phy_rm_sw <- check_rm_sw$tree_not_data
rm_dat_rm_sw <- check_rm_sw$data_not_tree
ctree_rm_sw <- drop.tip(tree, rm_phy_rm_sw)
#name.check(ctree_rm_sw, rm_sw)



SSX.rm_sw <- sum(round((rm_sw$rawFem - mean(rm_sw$rawFem))^2), 2)
s2.rm_sw <- var(rm_sw$rep_out_mass)
n.rm_sw <- length(rm_sw$rep_out_mass)
x.rm_sw <- seq(0.28, 3127.7, 1)
m.x.rm_sw <- mean(round(rm_sw$rawFem, 1))
se.rm_sw <- sqrt(s2.rm_sw*((1/n.rm_sw) + (((x.rm_sw - m.x.rm_sw)^2)/SSX.rm_sw)))
is.rm_sw <- qt(0.975, df = n.rm_sw - 2)
ii.rm_sw <- qt(0.025, df = n.rm_sw - 2)
ic.s.rm_sw <- se.rm_sw*is.rm_sw
ic.i.rm_sw <- se.rm_sw*ii.rm_sw
upper.i.rm_sw <- ((coef(mod_mass2)[1] + coef(mod_mass2)[3]) + (coef(mod_mass2)[2] + coef(mod_mass2)[5])*x.rm_sw) + ic.s.rm_sw
lower.i.rm_sw <- ((coef(mod_mass2)[1] + coef(mod_mass2)[3]) + (coef(mod_mass2)[2] + coef(mod_mass2)[5])*x.rm_sw) + ic.i.rm_sw



## mixed

rm_mx <- rawMass[rawMass$foraging.mode == "mixed", ]


check_rm_mx <- name.check(tree, rm_mx)
rm_phy_rm_mx <- check_rm_mx$tree_not_data
rm_dat_rm_mx <- check_rm_mx$data_not_tree
ctree_rm_mx <- drop.tip(tree, rm_phy_rm_mx)
#name.check(ctree_rm_mx, rm_mx)


SSX.rm_mx <- sum(round((rm_mx$rawFem - mean(rm_mx$rawFem))^2), 2)
s2.rm_mx <- var(rm_mx$rep_out_mass)
n.rm_mx <- length(rm_mx$rep_out_mass)
x.rm_mx <- seq(0.44, 222.1, 1)
m.x.rm_mx <- mean(round(rm_mx$rawFem, 1))
se.rm_mx <- sqrt(s2.rm_mx*((1/n.rm_mx) + (((x.rm_mx - m.x.rm_mx)^2)/SSX.rm_mx)))
is.rm_mx <- qt(0.975, df = n.rm_mx - 2)
ii.rm_mx <- qt(0.025, df = n.rm_mx - 2)
ic.s.rm_mx <- se.rm_mx*is.rm_mx
ic.i.rm_mx <- se.rm_mx*ii.rm_mx
upper.i.rm_mx.f <- ((coef(mod_mass2)[1] + coef(mod_mass2)[4]) + (coef(mod_mass2)[2] + coef(mod_mass2)[6])*x.rm_mx) + ic.s.rm_mx
lower.i.rm_mx.f <- ((coef(mod_mass2)[1] + coef(mod_mass2)[4]) + (coef(mod_mass2)[2] + coef(mod_mass2)[6])*x.rm_mx) + ic.i.rm_mx



par(mar = c(0, 5, 0, 0.5), mgp = c(3.2, 1, 0), las = 1)


plot(rawMass$rep_out_mass ~ rawMass$rawFem, type = "p", xlab = "mass (g)", ylab = "Hatchling mass x clutch size (g)", cex.lab = 1.5, pch = 21, bg = c("purple", "skyblue", "red")[as.numeric(rawMass$foraging.mode)], col = c("purple", "skyblue", "red")[as.numeric(rawMass$foraging.mode)])


polygon(c(rev(x.rm_mx), x.rm_mx), c(rev(lower.i.rm_mx.f), upper.i.rm_mx.f), border = FALSE, col = rgb(0.8, 0, 0, alpha = 0.3))
lines(x = x.rm_mx, y = ((coef(mod_mass2)[1] + coef(mod_mass2)[4]) + (coef(mod_mass2)[2] + coef(mod_mass2)[6])*x.rm_mx), col = "red", lwd = 2, lty = 15)

polygon(c(rev(x.rm_sw), x.rm_sw), c(rev(lower.i.rm_sw), upper.i.rm_sw), border = FALSE, col = rgb(0, 0, 1, alpha = 0.2))
lines(x = x.rm_sw, y = ((coef(mod_mass2)[1] + coef(mod_mass2)[3]) + (coef(mod_mass2)[2] + coef(mod_mass2)[5])*x.rm_sw), col = "skyblue", lwd = 2, lty = 2)

polygon(c(rev(x.rm_wf), x.rm_wf), c(rev(lower.i.rm_wf), upper.i.rm_wf), border = FALSE, col = rgb(0.8, 0, 0.7, alpha = 0.3))
lines(x = x.rm_wf, y = (coef(mod_mass2)[1] + coef(mod_mass2)[2]*x.rm_wf), col = "purple", lwd = 2, lty = 2)

mtext("(a)", side = 3, line = 1, at = -600, font = 2)
mtext("mass (g)", side = 1, line = 4, cex = 1)

legend("topleft", legend = c("mixed", "sit and wait", "widely foraging"), lty = c(15, 1, 2), bty = "n", pch = c(16, 16, 16), bg = c("black", "black", "black"), col = c("red", "skyblue", "purple"), cex = 1)





## Scaled mass index

## widely foraging

wf <- dat_full[dat_full$foraging.mode == "widely foraging", ]

check_wf <- name.check(tree, wf)
rm_phy_wf <- check_wf$tree_not_data
rm_dat_wf <- check_wf$data_not_tree
ctree_wf <- drop.tip(tree, rm_phy_wf)
#name.check(ctree_wf, wf)


SSX.wf <- sum(round((wf$M_females - mean(wf$M_females))^2), 2)
s2.wf <- var(wf$rep_out)
n.wf <- length(wf$rep_out)
x.wf <- seq(1.2, 31.8, 1)
m.x.wf <- mean(round(wf$M_females,1))
se.wf <- sqrt(s2.wf*((1/n.wf) + (((x.wf - m.x.wf)^2)/SSX.wf)))
is.wf <- qt(0.975, df = n.wf - 2)
ii.wf <- qt(0.025, df = n.wf - 2)
ic.s.wf <- se.wf*is.wf
ic.i.wf <- se.wf*ii.wf
upper.i <- (coef(full_mod2)[1] + coef(full_mod2)[2]*x.wf) + ic.s.wf
lower.i <- (coef(full_mod2)[1] + coef(full_mod2)[2]*x.wf) + ic.i.wf



## sit and wait


sw <- dat_full[dat_full$foraging.mode == "sit and wait", ]

check_sw <- name.check(tree, sw)
rm_phy_sw <- check_sw$tree_not_data
rm_dat_sw <- check_sw$data_not_tree
ctree_sw <- drop.tip(tree, rm_phy_sw)
#name.check(ctree_sw, sw)



SSX.sw <- sum(round((sw$M_females - mean(sw$M_females))^2), 2)
s2.sw <- var(sw$rep_out)
n.sw <- length(sw$rep_out)
x.sw <- seq(1.3, 30.8, 1)
m.x.sw <- mean(round(sw$M_females, 1))
se.sw <- sqrt(s2.sw*((1/n.sw) + (((x.sw - m.x.sw)^2)/SSX.sw)))
is.sw <- qt(0.975, df = n.sw - 2)
ii.sw <- qt(0.025, df = n.sw - 2)
ic.s.sw <- se.sw*is.sw
ic.i.sw <- se.sw*ii.sw
upper.i.sw <- ((coef(full_mod2)[1] + coef(full_mod2)[3]) + (coef(full_mod2)[2] + coef(full_mod2)[5])*x.sw) + ic.s.sw
lower.i.sw <- ((coef(full_mod2)[1] + coef(full_mod2)[3]) + (coef(full_mod2)[2] + coef(full_mod2)[5])*x.sw) + ic.i.sw


## mixed

mx <- dat_full[dat_full$foraging.mode == "mixed", ]


check_mx <- name.check(tree, mx)
rm_phy_mx <- check_mx$tree_not_data
rm_dat_mx <- check_mx$data_not_tree
ctree_mx <- drop.tip(tree, rm_phy_mx)
#name.check(ctree_mx, mx)


SSX.mx <- sum(round((mx$M_females - mean(mx$M_females))^2), 2)
s2.mx <- var(mx$rep_out)
n.mx <- length(mx$rep_out)
x.mx <- seq(1.0, 27.7, 1)
m.x.mx <- mean(round(mx$M_females, 1))
se.mx <- sqrt(s2.mx*((1/n.mx) + (((x.mx - m.x.mx)^2)/SSX.mx)))
is.mx <- qt(0.975, df = n.mx - 2)
ii.mx <- qt(0.025, df = n.mx - 2)
ic.s.mx <- se.mx*is.mx
ic.i.mx <- se.mx*ii.mx
upper.i.mx.f <- ((coef(full_mod2)[1] + coef(full_mod2)[4]) + (coef(full_mod2)[2] + coef(full_mod2)[6])*x.mx) + ic.s.mx
lower.i.mx.f <- ((coef(full_mod2)[1] + coef(full_mod2)[4]) + (coef(full_mod2)[2] + coef(full_mod2)[6])*x.mx) + ic.i.mx


par(mar = c(0, 4, 0, 0.5), mgp = c(2.8, 1, 0), las = 1)


plot(dat_full$rep_out ~ dat_full$M_females, type = "p", xlab = " ", ylab = "Hatchling mass x clutch size (g)", cex.lab = 1.5, pch = 21, bg = c("purple", "skyblue", "red")[as.numeric(dat_full$foraging.mode)], col = c("purple", "skyblue", "red")[as.numeric(dat_full$foraging.mode)])

polygon(c(rev(x.mx), x.mx), c(rev(lower.i.mx.f), upper.i.mx.f), border = FALSE, col = rgb(0.8, 0, 0, alpha = 0.3))
lines(x = x.mx, y = ((coef(full_mod2)[1] + coef(full_mod2)[4]) + (coef(full_mod2)[2] + coef(full_mod2)[6])*x.mx), col = "red", lwd = 2, lty = 15)

polygon(c(rev(x.sw), x.sw), c(rev(lower.i.sw), upper.i.sw), border = FALSE, col = rgb(0, 0, 1, alpha = 0.2))
lines(x = x.sw, y = ((coef(full_mod2)[1] + coef(full_mod2)[3]) + (coef(full_mod2)[2] + coef(full_mod2)[5])*x.sw), col = "skyblue", lwd = 2, lty = 2)

polygon(c(rev(x.wf), x.wf), c(rev(lower.i), upper.i), border = FALSE, col = rgb(0.8, 0, 0.7, alpha = 0.3))
lines(x = x.wf, y = coef(full_mod2)[1] + coef(full_mod2)[2]*x.wf, col = "purple", lwd = 2, lty = 2)

mtext("(b)", side = 3, line = 1, at = -2, font = 2)
mtext(expression(hat(M) (g)), side = 1, line = 4, cex = 1)

legend("topleft", legend = c("mixed", "sit-and-wait foraging", "widely foraging"), lty = c(15, 1, 2), bty = "n", pch = c(16, 16, 16), bg = c("black", "black", "black"), col = c("red", "skyblue", "purple"), cex = 1)

Figure 3. Effects of maternal mass and foraging mode on the evolution of reproductive output of lizards, as determined by phylogenetic generalized least squares analysis. (a) Difference in reproductive output among foraging modes assuming body mass as a proxy for energy. (b) Difference in reproductive output among foraging modes assuming the scaled mass index as a proxy for energy.




Ancenstral state reconstruction

dev.off()
  null device 
            1
set.seed(1094)


## Ancestral character state

anc_st <- cdat[ , c(25, 23)]
anc_st[486, ] <- c("Sphenodon_punctatus", "sit and wait")
rownames(anc_st)[486] <- "Sphenodon_punctatus"

check_st <- name.check(tree, anc_st)
rm_phy_st <- check_st$tree_not_data
rm_dat_st <- check_st$data_not_tree
ctree_st <- drop.tip(tree, rm_phy_st)

anc_st <- subset(anc_st, subset = rownames(anc_st) %in% ctree_st$tip, select = names(anc_st)[1:2])
#name.check(ctree, anc_st)

forg <- as.matrix(anc_st)[,2]
head(forg)
       Acanthocercus_atricollis Acanthodactylus_beershebensis 
                 "sit and wait"             "widely foraging" 
      Acanthodactylus_boskianus    Acanthodactylus_erythrurus 
              "widely foraging"                       "mixed" 
       Acanthodactylus_pardalis    Acanthodactylus_schreiberi 
                 "sit and wait"             "widely foraging"
## Reproductive output

anc_st_rp <- cdat[ , c(29, 23)]
anc_st_rp[486, 1] <- 0
anc_st_rp[486, 2] <- "sit and wait"
rownames(anc_st_rp)[486] <- "Sphenodon_punctatus"

check_rp <- name.check(tree, anc_st_rp)
rm_phy_rp <- check_rp$tree_not_data
rm_dat_rp <- check_rp$data_not_tree
ctree_rp <- drop.tip(tree, rm_phy_rp)

anc_st_rp <- subset(anc_st_rp, subset = rownames(anc_st_rp) %in% ctree_rp$tip, select = names(anc_st_rp)[1:2])
#name.check(ctree, anc_st_rp)
anc_st_rp$rep_out <- as.numeric(anc_st_rp$rep_out)
anc_st_rp$foraging.mode <- as.factor(anc_st_rp$foraging.mode)

anc_st_rp_bars <- anc_st_rp[1]
anc_st_rp_bars <- as.matrix(anc_st_rp_bars)[ , 1]
head(anc_st_rp_bars)
       Acanthocercus_atricollis Acanthodactylus_beershebensis 
                      11.268193                      3.720936 
      Acanthodactylus_boskianus    Acanthodactylus_erythrurus 
                       3.284710                      2.879015 
       Acanthodactylus_pardalis    Acanthodactylus_schreiberi 
                       2.845583                      2.252075
## Reproductive output for widely foraging and sit-and-wait species

anc_st_rp2 <- anc_st_rp[anc_st_rp$foraging.mode != "mixed", ]

check_rp2 <- name.check(tree, anc_st_rp2)
rm_phy_rp2 <- check_rp2$tree_not_data
rm_dat_rp2 <- check_rp2$data_not_tree
ctree_rp2 <- drop.tip(tree, rm_phy_rp2)

anc_st_rp2 <- subset(anc_st_rp2, subset = rownames(anc_st_rp2) %in% ctree_rp2$tip, select = names(anc_st_rp2)[1:2])
#name.check(ctree_rp2, anc_st_rp2)
anc_st_rp2$rep_out <- as.numeric(anc_st_rp2$rep_out)
anc_st_rp2$foraging.mode <- as.factor(anc_st_rp2$foraging.mode)

anc_st_rp2_bars <- anc_st_rp2[1]
anc_st_rp2_bars <- as.matrix(anc_st_rp2_bars)[ , 1]
head(anc_st_rp2_bars)
       Acanthocercus_atricollis Acanthodactylus_beershebensis 
                     11.2681934                     3.7209360 
      Acanthodactylus_boskianus      Acanthodactylus_pardalis 
                      3.2847103                     2.8455830 
     Acanthodactylus_schreiberi          Acontias_gariepensis 
                      2.2520747                     0.7391738
## Use of function ace (for comparison)


fitER <- ace(forg, ctree_st, model = "ER", type = "discrete")
fitARD <- ace(forg, ctree_st, model = "ARD", type = "discrete")
fitSR <- ace(forg, ctree_st, model = "SYM", type = "discrete")


fitER
  
      Ancestral Character Estimation
  
  Call: ace(x = forg, phy = ctree_st, type = "discrete", model = "ER")
  
      Log-likelihood: -335.5832 
  
  Rate index matrix:
                  mixed sit and wait widely foraging
  mixed               .            1               1
  sit and wait        1            .               1
  widely foraging     1            1               .
  
  Parameter estimates:
   rate index estimate std-err
            1   0.0032   3e-04
  
  Scaled likelihoods at the root (type '...$lik.anc' to get them for all nodes):
            mixed    sit and wait widely foraging 
        0.1631016       0.4946075       0.3422909
fitARD
  
      Ancestral Character Estimation
  
  Call: ace(x = forg, phy = ctree_st, type = "discrete", model = "ARD")
  
      Log-likelihood: -320.1167 
  
  Rate index matrix:
                  mixed sit and wait widely foraging
  mixed               .            3               5
  sit and wait        1            .               6
  widely foraging     2            4               .
  
  Parameter estimates:
   rate index estimate std-err
            1   0.0062  0.0011
            2   0.0045  0.0010
            3   0.0133  0.0040
            4   0.0022  0.0007
            5   0.0205  0.0051
            6   0.0000  0.0007
  
  Scaled likelihoods at the root (type '...$lik.anc' to get them for all nodes):
            mixed    sit and wait widely foraging 
        0.3347980       0.3811181       0.2840839
fitSR
  
      Ancestral Character Estimation
  
  Call: ace(x = forg, phy = ctree_st, type = "discrete", model = "SYM")
  
      Log-likelihood: -333.6765 
  
  Rate index matrix:
                  mixed sit and wait widely foraging
  mixed               .            1               2
  sit and wait        1            .               3
  widely foraging     2            3               .
  
  Parameter estimates:
   rate index estimate std-err
            1   0.0041   6e-04
            2   0.0041   7e-04
            3   0.0023   4e-04
  
  Scaled likelihoods at the root (type '...$lik.anc' to get them for all nodes):
            mixed    sit and wait widely foraging 
        0.2145978       0.4773200       0.3080822
head(round(fitER$lik.anc, 3))
      mixed sit and wait widely foraging
  487 0.163        0.495           0.342
  488 0.018        0.582           0.399
  489 0.005        0.954           0.041
  490 0.010        0.969           0.021
  491 0.017        0.953           0.030
  492 0.069        0.921           0.011
head(round(fitARD$lik.anc, 3))
      mixed sit and wait widely foraging
  487 0.335        0.381           0.284
  488 0.283        0.362           0.356
  489 0.074        0.895           0.031
  490 0.093        0.893           0.013
  491 0.132        0.853           0.015
  492 0.138        0.855           0.007
head(round(fitSR$lik.anc, 3))
      mixed sit and wait widely foraging
  487 0.215        0.477           0.308
  488 0.060        0.566           0.374
  489 0.017        0.953           0.030
  490 0.024        0.962           0.014
  491 0.036        0.943           0.021
  492 0.089        0.904           0.007




Table S2.

ic <- data.frame(Model = c("ER", "ARD", "SYM"), AIC = c(AIC(fitER), AIC(fitARD), AIC(fitSR)), stringsAsFactors = FALSE)

kable(ic, digits = 3, row.names = FALSE) %>%
  kable_styling(bootstrap_options = "striped", full_width = FALSE, position = "center")
Model AIC
ER 673.166
ARD 652.233
SYM 673.353




## Overlaying posterior probabilities on the tree. This is done for the best fit model.


colors <- c("red", "skyblue", "purple")
cols <- setNames(colors[1:length(unique(forg))], sort(unique(forg)))


plotTree(ctree_st, type = "fan", fsize = 0.2, ftype = "i", part = 0.98)
t <- max(nodeHeights(ctree_st))
tick.spacing <- 25
min.tick <- 0
scale <- axis(1, pos = -5, at = seq(t, min.tick, by = -tick.spacing), cex.axis = 0.5, labels = FALSE)

for(i in 1:length(scale)){
    a1 <- 0
    a2 <- 2*pi
    draw.arc(0, 0, radius = scale[i], a1, a2, lwd = 1,
        col = make.transparent("blue", 0.15))
}

text(115, 45, "1", cex = 1)
text(32, 90, "2", cex = 1)
text(-105, 25, "3", cex = 1)
text(-135, -70, "4", cex = 1)
text(-45, -65, "5", cex = 1)
text(-8, -95, "6", cex = 1)
text(scale, rep(-23, length(scale)), t - scale, cex = 0.6)
text(mean(scale), -38, "time (mya)", cex = 0.7)

nodelabels(node = 1:ctree_st$Nnode + Ntip(ctree_st), pie = fitARD$lik.anc, piecol = cols, cex = 0.2)
tiplabels(pie = to.matrix(forg[ctree_st$tip.label], unique(sort(forg))), piecol = cols, cex = 0.1)
add.simmap.legend(colors = cols, prompt = FALSE, x= 0.9*par()$usr[1], y = -max(nodeHeights(ctree_st)), fsize = 0.8)

ASER <- make.simmap(ctree_st, forg,  model = "ER", nsim = 1, pi = "estimated", Q = "mcmc")
  Running MCMC burn-in. Please wait....
  Running 100 generations of MCMC, sampling every 100 generations.
  Please wait....
  
  make.simmap is simulating with a sample of Q from
  the posterior distribution
  
  Mean Q from the posterior is
  Q =
                         mixed sit and wait widely foraging
  mixed           -0.006716348  0.003358174     0.003358174
  sit and wait     0.003358174 -0.006716348     0.003358174
  widely foraging  0.003358174  0.003358174    -0.006716348
  and (mean) root node prior probabilities
  pi =
            mixed    sit and wait widely foraging 
        0.3333333       0.3333333       0.3333333
ASARD <- make.simmap(ctree_st, forg,  model = "ARD", nsim = 1, pi = "estimated", Q = "mcmc")
  Running MCMC burn-in. Please wait....
  Running 100 generations of MCMC, sampling every 100 generations.
  Please wait....
  
  make.simmap is simulating with a sample of Q from
  the posterior distribution
  
  Mean Q from the posterior is
  Q =
                       mixed sit and wait widely foraging
  mixed           -4.2957510    3.8697964       0.4259546
  sit and wait     0.6786649   -2.1882412       1.5095763
  widely foraging  0.5975036    0.9830754      -1.5805790
  and (mean) root node prior probabilities
  pi =
            mixed    sit and wait widely foraging 
        0.1291999       0.4275951       0.4432051
ASSR <- make.simmap(ctree_st, forg,  model = "SYM", nsim = 1, pi = "estimated", Q = "mcmc")
  Running MCMC burn-in. Please wait....
  Running 100 generations of MCMC, sampling every 100 generations.
  Please wait....
  
  make.simmap is simulating with a sample of Q from
  the posterior distribution
  
  Mean Q from the posterior is
  Q =
                        mixed sit and wait widely foraging
  mixed           -1.56999075   0.01032159       1.5596692
  sit and wait     0.01032159  -0.02980089       0.0194793
  widely foraging  1.55966916   0.01947930      -1.5791485
  and (mean) root node prior probabilities
  pi =
            mixed    sit and wait widely foraging 
        0.3333333       0.3333333       0.3333333




Table S3.

as <- data.frame(Model = c("modER", "modARD", "modSYM"), AIC = c(AIC(ASER), AIC(ASARD), AIC(ASSR)), stringsAsFactors = FALSE)

kable(as, digits = 3, row.names = FALSE) %>%
  kable_styling(bootstrap_options = "striped", full_width = FALSE, position = "center")
Model AIC
modER 685.665
modARD 969.306
modSYM 947.630




## Simulate 1000 stochastic character maps using empirical Bayes method


ASER <- make.simmap(ctree_st, forg,  model = "ER", nsim = 1000, pi = "estimated", Q = "mcmc")
  Running MCMC burn-in. Please wait....
  Running 1e+05 generations of MCMC, sampling every 100 generations.
  Please wait....
  
  make.simmap is simulating with a sample of Q from
  the posterior distribution
  
  Mean Q from the posterior is
  Q =
                         mixed sit and wait widely foraging
  mixed           -0.006463354  0.003231677     0.003231677
  sit and wait     0.003231677 -0.006463354     0.003231677
  widely foraging  0.003231677  0.003231677    -0.006463354
  and (mean) root node prior probabilities
  pi =
            mixed    sit and wait widely foraging 
        0.3333333       0.3333333       0.3333333
pd <- summary(ASER, plot = FALSE)
pd
  1000 trees with a mapped discrete character with states:
   mixed, sit and wait, widely foraging 
  
  trees have 126.77 changes between states on average
  
  changes are of the following types:
       mixed,sit and wait mixed,widely foraging sit and wait,mixed
  x->y               8.05                  8.06             34.261
       sit and wait,widely foraging widely foraging,mixed
  x->y                       24.335                25.801
       widely foraging,sit and wait
  x->y                       26.263
  
  mean total time spent in each state is:
              mixed sit and wait widely foraging    total
  raw  1.648428e+03 1.033519e+04    7900.6573103 19884.28
  prop 8.290106e-02 5.197671e-01       0.3973319     1.00
save <- countSimmap(ASER)


min(save$Tr[,1])
  [1] 105
write.table(pd$ace, "likelihood_nodes_ASER.csv", sep = ",")



## What happens if the tuatara is excluded?


forgt <- forg[-486]


ASERT <- make.simmap(ctree, forgt,  model = "ER", nsim = 1000, pi = "estimated", Q = "mcmc")
  Running MCMC burn-in. Please wait....
  Running 1e+05 generations of MCMC, sampling every 100 generations.
  Please wait....
  
  make.simmap is simulating with a sample of Q from
  the posterior distribution
  
  Mean Q from the posterior is
  Q =
                         mixed sit and wait widely foraging
  mixed           -0.006425791  0.003212896     0.003212896
  sit and wait     0.003212896 -0.006425791     0.003212896
  widely foraging  0.003212896  0.003212896    -0.006425791
  and (mean) root node prior probabilities
  pi =
            mixed    sit and wait widely foraging 
        0.3333333       0.3333333       0.3333333
pdt <- summary(ASERT, plot = FALSE)
pdt
  1000 trees with a mapped discrete character with states:
   mixed, sit and wait, widely foraging 
  
  trees have 124.251 changes between states on average
  
  changes are of the following types:
       mixed,sit and wait mixed,widely foraging sit and wait,mixed
  x->y              7.342                 7.678             33.895
       sit and wait,widely foraging widely foraging,mixed
  x->y                       24.041                25.419
       widely foraging,sit and wait
  x->y                       25.876
  
  mean total time spent in each state is:
              mixed sit and wait widely foraging    total
  raw  1.588090e+03 1.012479e+04     7818.285105 19531.16
  prop 8.131054e-02 5.183915e-01        0.400298     1.00
savet <- countSimmap(ASERT)


min(savet$Tr[,1])
  [1] 102
write.table(pdt$ace, "likelihood_nodes_ASERT.csv", sep = ",")


## Phylogenetic signal of foraging modes


anc_st_rp2$Binomial <- rownames(anc_st_rp2)

phylo.d(anc_st_rp2, ctree_rp2, Binomial, binvar = foraging.mode, permut = 5000)
  
  Calculation of D statistic for the phylogenetic structure of a binary variable
  
    Data :  data
    Binary variable :  foraging.mode
    Counts of states:  widely foraging = 202
                       sit and wait = 223
    Phylogeny :  phy
    Number of permutations :  5000
  
  Estimated D :  -0.144269
  Probability of E(D) resulting from no (random) phylogenetic structure :  0
  Probability of E(D) resulting from Brownian phylogenetic structure    :  0.8004
## Estimating ancestral character state only for widely foraging and sit-and-wait species


forg2 <- as.matrix(anc_st_rp2)[ , 2]
head(forg2)
       Acanthocercus_atricollis Acanthodactylus_beershebensis 
                 "sit and wait"             "widely foraging" 
      Acanthodactylus_boskianus      Acanthodactylus_pardalis 
              "widely foraging"                "sit and wait" 
     Acanthodactylus_schreiberi          Acontias_gariepensis 
              "widely foraging"             "widely foraging"
ASER2 <- make.simmap(ctree_rp2, forg2, model = "ER", nsim = 1, pi = "estimated", Q = "mcmc")
  Running MCMC burn-in. Please wait....
  Running 100 generations of MCMC, sampling every 100 generations.
  Please wait....
  
  make.simmap is simulating with a sample of Q from
  the posterior distribution
  
  Mean Q from the posterior is
  Q =
                  sit and wait widely foraging
  sit and wait    -0.002167574     0.002167574
  widely foraging  0.002167574    -0.002167574
  and (mean) root node prior probabilities
  pi =
     sit and wait widely foraging 
              0.5             0.5
ASARD2 <- make.simmap(ctree_rp2, forg2, model = "ARD", nsim = 1, pi = "estimated", Q = "mcmc")
  Running MCMC burn-in. Please wait....
  Running 100 generations of MCMC, sampling every 100 generations.
  Please wait....
  
  make.simmap is simulating with a sample of Q from
  the posterior distribution
  
  Mean Q from the posterior is
  Q =
                  sit and wait widely foraging
  sit and wait    -0.002013939     0.002013939
  widely foraging  0.007239599    -0.007239599
  and (mean) root node prior probabilities
  pi =
     sit and wait widely foraging 
        0.7823601       0.2176399
ASSR2 <- make.simmap(ctree_rp2, forg2, model = "SYM", nsim = 1, pi = "estimated", Q = "mcmc")
  Running MCMC burn-in. Please wait....
  Running 100 generations of MCMC, sampling every 100 generations.
  Please wait....
  
  make.simmap is simulating with a sample of Q from
  the posterior distribution
  
  Mean Q from the posterior is
  Q =
                  sit and wait widely foraging
  sit and wait     -0.00306315      0.00306315
  widely foraging   0.00306315     -0.00306315
  and (mean) root node prior probabilities
  pi =
     sit and wait widely foraging 
              0.5             0.5




Table S4.

as2 <- data.frame(Model = c("modER2", "modARD2", "modSR2"), AIC = c(AIC(ASER2), AIC(ASARD2), AIC(ASSR2)), stringsAsFactors = FALSE)
kable(as2, digits = 3, row.names = FALSE) %>%
  kable_styling(bootstrap_options = "striped", full_width = FALSE, position = "center")
Model AIC
modER2 297.869
modARD2 309.341
modSR2 295.076



## simulate 1000 stochastic character maps using empirical Bayes method (excluding mixed foragers)


ASSYM2 <- make.simmap(ctree_rp2, forg2, model = "SYM", nsim = 1000, pi = "estimated", Q = "mcmc")
  Running MCMC burn-in. Please wait....
  Running 1e+05 generations of MCMC, sampling every 100 generations.
  Please wait....
  
  make.simmap is simulating with a sample of Q from
  the posterior distribution
  
  Mean Q from the posterior is
  Q =
                  sit and wait widely foraging
  sit and wait    -0.003025366     0.003025366
  widely foraging  0.003025366    -0.003025366
  and (mean) root node prior probabilities
  pi =
     sit and wait widely foraging 
              0.5             0.5
pd2 <- summary(ASSYM2, plot = FALSE)
pd2
  1000 trees with a mapped discrete character with states:
   sit and wait, widely foraging 
  
  trees have 52.577 changes between states on average
  
  changes are of the following types:
       sit and wait,widely foraging widely foraging,sit and wait
  x->y                       26.137                        26.44
  
  mean total time spent in each state is:
       sit and wait widely foraging    total
  raw  1.010405e+04    7712.4930593 17816.54
  prop 5.671161e-01       0.4328839     1.00
save2 <- countSimmap(ASSYM2)

min(save2$Tr[,1])
  [1] 42
write.table(pd2$ace, "likelihood_nodes_ASSYM2.csv", sep = ",")



## simulate 1000 stochastic character maps using empirical Bayes method (excluding mixed foragers and the tuatara)


forg2T <- cdat[ , c(25, 23)]
forg2T <- forg2T[forg2T$foraging.mode != "mixed", ]

check_st2 <- name.check(tree, forg2T)
rm_phy_st2 <- check_st2$tree_not_data
rm_dat_st2 <- check_st2$data_not_tree
ctree_st2 <- drop.tip(tree, rm_phy_st2)

forg2T <- subset(forg2T, subset = rownames(forg2T) %in% ctree_st2$tip, select = names(forg2T)[1:2])
#name.check(ctree_st2, forg2T)

forg2T <- as.matrix(forg2T)[,2]
head(forg2T)
       Acanthocercus_atricollis Acanthodactylus_beershebensis 
                 "sit and wait"             "widely foraging" 
      Acanthodactylus_boskianus      Acanthodactylus_pardalis 
              "widely foraging"                "sit and wait" 
     Acanthodactylus_schreiberi          Acontias_gariepensis 
              "widely foraging"             "widely foraging"
ASSYMT <- make.simmap(ctree_st2, forg2T, model = "SYM", nsim = 1000, pi = "estimated", Q = "mcmc")
  Running MCMC burn-in. Please wait....
  Running 1e+05 generations of MCMC, sampling every 100 generations.
  Please wait....
  
  make.simmap is simulating with a sample of Q from
  the posterior distribution
  
  Mean Q from the posterior is
  Q =
                  sit and wait widely foraging
  sit and wait    -0.002965521     0.002965521
  widely foraging  0.002965521    -0.002965521
  and (mean) root node prior probabilities
  pi =
     sit and wait widely foraging 
              0.5             0.5
pd2T <- summary(ASSYMT, plot = FALSE)
pd2T
  1000 trees with a mapped discrete character with states:
   sit and wait, widely foraging 
  
  trees have 51.071 changes between states on average
  
  changes are of the following types:
       sit and wait,widely foraging widely foraging,sit and wait
  x->y                       24.979                       26.092
  
  mean total time spent in each state is:
       sit and wait widely foraging    total
  raw   9788.337854     7675.089625 17463.43
  prop     0.560505        0.439495     1.00
save2T <- countSimmap(ASSYMT)

min(save2T$Tr[,1])
  [1] 42
write.table(pd2T$ace, "likelihood_nodes_ASSYMT.csv", sep = ",")


## Plotting trees

colors <- c("red", "skyblue", "purple")
cols <- setNames(colors[1:length(unique(forg))], sort(unique(forg)))
colorbars <-setNames(c("red", "skyblue", "purple"), c("mixed", "sit and wait", "widely foraging"))
ss <-getStates(ASER, "tips")
barcols <- setNames(sapply(ss, function(x, y) y[which(names(y) == x)], y = colorbars), names(ss))

i <- 1:1000

plotTree.wBars(ASER[[sample(i, 1)]], log1p(anc_st_rp_bars), type = "fan", method = "plotSimmap", colors = cols, scale = 4, fsize = 0.2, tip.labels = FALSE, part = 0.98, col = barcols, border = FALSE)
t1 <- max(nodeHeights(ctree_st))
tick.spacing <- 25
min.tick <- 0
scale <- axis(1, pos = -5, at = seq(t1, min.tick, by = -tick.spacing), cex.axis = 0.5, labels = FALSE)

for(i in 1:length(scale)){
    a1 <- 0
    a2 <- 2*pi
    draw.arc(0, 0, radius = scale[i], a1, a2, lwd = 1,
        col = make.transparent("blue", 0.15))
}

text(119, 45, "1", cex = 1)
text(32, 90, "2", cex = 1)
text(-105, 25, "3", cex = 1)
text(-135, -70, "4", cex = 1)
text(-45, -65, "5", cex = 1)
text(-8, -95, "6", cex = 1)
text(scale, rep(-23, length(scale)), t1 - scale, cex = 0.9)
text(mean(scale), -40, "time (mya)", cex = 0.9)
nodelabels(pie = pd$ace, piecol = cols, cex = 0.26)
add.simmap.legend(colors = cols, prompt = FALSE, x = -120, y = 140, fsize = 0.9)

Figure 2. Random sample of stochastic character maps depicting the evolution of foraging mode in 485 species of lizards. Bars at the tips of the phylogeny represent log-transformed values of reproductive output for all lizards, but not the outgroup,Sphenodon punctatus. Pie charts on internal nodes represent posterior probability estimates resulting from 1,000 simulations of the character histories. Major clades are enumerated as follows: 1) Gekkota, 2) Scincoidea, 3) Lacertoidea, 4) Anguimorpha, 5) Toxicofera, and 6) Iguania. Lizard photos by Mark O’Shea.

## Ancestral state reconstruction excluding the tuatara

colst <- setNames(colors[1:length(unique(forgt))], sort(unique(forgt)))
sst <-getStates(ASERT, "tips")
colorbars <-setNames(c("red", "skyblue", "purple"), c("mixed", "sit and wait", "widely foraging"))
barcolst <- setNames(sapply(sst, function(x, y) y[which(names(y) == x)], y = colorbars), names(sst))

iT <- 1:1000

plotTree.wBars(ASERT[[sample(iT, 1)]], log1p(anc_st_rp_bars[-486]), type = "fan", method = "plotSimmap", colors = colst, scale = 4, fsize = 0.2, tip.labels = FALSE, part = 0.98, col = barcolst, border = FALSE)
t2 <- max(nodeHeights(ctree))
tick.spacing <- 25
min.tick <- 0
scale <- axis(1, pos = -5, at = seq(t2, min.tick, by = -tick.spacing), cex.axis = 0.5, labels = FALSE)

for(i in 1:length(scale)){
    a1 <- 0
    a2 <- 2*pi
    draw.arc(0, 0, radius = scale[i], a1, a2, lwd = 1,
        col = make.transparent("blue", 0.15))
}

text(55, 20, "1", cex = 1)
text(10, 20, "2", cex = 1)
text(-30, 8, "3", cex = 1)
text(-72, -38, "4", cex = 1)
text(-11, -24, "5", cex = 1)
text(-8, -40, "6", cex = 1)
text(scale, rep(-20, length(scale)), t2 - scale, cex = 0.7)
text(mean(scale), -35, "time (mya)", cex = 0.9)
nodelabels(pie = pdt$ace, piecol = colst, cex = 0.3)
add.simmap.legend(colors = colst, prompt = FALSE, x = -80, y = 90, fsize = 0.9)




## Internal nodes numbering

plotTree.wBars(ASER[[sample(i, 1)]], log1p(anc_st_rp_bars), type = "fan", method = "plotSimmap", colors = cols, scale = 4, fsize = 0.2, tip.labels = FALSE, part = 0.98, col = barcols, border = FALSE)
t1 <- max(nodeHeights(ctree_st))
tick.spacing <- 25
min.tick <- 0
scale <- axis(1, pos = -5, at = seq(t1, min.tick, by = -tick.spacing), cex.axis = 0.5, labels = FALSE)

for(i in 1:length(scale)){
    a1 <- 0
    a2 <- 2*pi
    draw.arc(0, 0, radius = scale[i], a1, a2, lwd = 1,
        col = make.transparent("blue", 0.15))
}

text(115, 45, "1", cex = 1)
text(32, 90, "2", cex = 1)
text(-105, 25, "3", cex = 1)
text(-135, -70, "4", cex = 1)
text(-45, -65, "5", cex = 1)
text(-8, -95, "6", cex = 1)
text(scale, rep(-23, length(scale)), t1 - scale, cex = 0.6)
text(mean(scale), -40, "time (mya)", cex = 0.7)
nodelabels(bg = "white", cex = 0.5)

## Standard tree

plot(ASER[[sample(i, 1)]], colors = cols, ftype = "off")

## Standard tree

plot(ASER[[sample(i, 1)]], colors = cols, ftype = "off")
nodelabels(bg = "white", cex = 0.5)

## Ancestral character state reconstruction without mixed foraging species


colors2 <- c("skyblue", "purple")
cols2 <- setNames(colors2[1:length(unique(forg2))], sort(unique(forg2)))
colorbars2 <-setNames(c("skyblue", "purple"), c("sit and wait", "widely foraging"))
ss2 <- getStates(ASSYM2, "tips")
barcols2 <- setNames(sapply(ss2, function(x, y) y[which(names(y) == x)], y = colorbars2), names(ss2))

i2 <- 1:1000

plotTree.wBars(ASSYM2[[sample(i2, 1)]], log1p(anc_st_rp2_bars), type = "fan", method = "plotSimmap", colors = cols2, scale = 4, fsize = 0.2, tip.labels = FALSE, part = 0.98, col = barcols2, border = FALSE)
t3 <- max(nodeHeights(ctree_rp2))
tick.spacing <- 25
min.tick <- 0
scale <- axis(1, pos = -5, at = seq(t3, min.tick, by = -tick.spacing), cex.axis = 0.5, labels = FALSE)

for(i in 1:length(scale)){
    a1 <- 0
    a2 <- 2*pi
    draw.arc(0, 0, radius = scale[i], a1, a2, lwd = 1,
        col = make.transparent("blue", 0.15))
}

text(115, 45, "1", cex = 1)
text(45, 85, "2", cex = 1)
text(-105, 35, "3", cex = 1)
text(-137, -68, "4", cex = 1)
text(-45, -65, "5", cex = 1)
text(-8, -95, "6", cex = 1)
text(scale, rep(-23, length(scale)), t3 - scale, cex = 0.9)
text(mean(scale), -40, "time (mya)", cex = 0.9)
nodelabels(pie = pd2$ace, piecol = cols2, cex = 0.26)
add.simmap.legend(colors = cols2, prompt = FALSE, x = -120, y = 140, fsize = 0.8)

## Excluding the tuatara

i2T <- 1:1000

colors2T <- c("skyblue", "purple")
cols2T <- setNames(colors2T[1:length(unique(forg2T))], sort(unique(forg2T)))
colorbars2T <-setNames(c("skyblue", "purple"), c("sit and wait", "widely foraging"))
ss2T <- getStates(ASSYMT, "tips")
barcols2T <- setNames(sapply(ss2T, function(x, y) y[which(names(y) == x)], y = colorbars2T), names(ss2T))


plotTree.wBars(ASSYMT[[sample(i2T, 1)]], log1p(anc_st_rp2_bars[-425]), type = "fan", method = "plotSimmap", colors = cols2T, scale = 4, fsize = 0.2, tip.labels = FALSE, part = 0.98, col = barcols2T, border = FALSE)
t4 <- max(nodeHeights(ctree_st2))
tick.spacing <- 25
min.tick <- 0
scale <- axis(1, pos = -5, at = seq(t4, min.tick, by = -tick.spacing), cex.axis = 0.5, labels = FALSE)

for(i in 1:length(scale)){
    a1 <- 0
    a2 <- 2*pi
    draw.arc(0, 0, radius = scale[i], a1, a2, lwd = 1,
        col = make.transparent("blue", 0.15))
}

text(55, 20, "1", cex = 1)
text(14, 20, "2", cex = 1)
text(-30, 8, "3", cex = 1)
text(-72, -35, "4", cex = 1)
text(-11, -24, "5", cex = 1)
text(-8, -40, "6", cex = 1)
text(scale, rep(-20, length(scale)), t4 - scale, cex = 0.7)
text(mean(scale), -40, "time (mya)", cex = 0.9)
nodelabels(pie = pd2T$ace, piecol = cols2T, cex = 0.26)
add.simmap.legend(colors = cols2T, prompt = FALSE, x = -70, y = 80, fsize = 0.8)